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Both  circadian  rhythmicity  and  sleep  play  significant  roles  in  the 
regulation  of  plasma  cortisol  concentration  by  the  hypothalamo- 
pituitary-adrenal  (HPA)  axis.  Numerous  studies  have  found  links 
between  sleep  and  changes  in  cortisol  concentration,  but  the  implica¬ 
tions  of  these  results  have  remained  largely  qualitative.  In  this  article, 
we  present  a  quantitative  phenomenological  model  to  describe  the 
effects  of  different  sleep  durations  on  cortisol  concentration.  We 
constructed  the  proposed  model  by  incorporating  the  circadian  and 
sleep  allostatic  effects  on  cortisol  concentration,  the  pulsatile  nature  of 
cortisol  secretion,  and  cortisol’s  negative  autoregulation  of  its  own 
production  and  validated  its  performance  on  three  study  groups  that 
experienced  four  distinct  sleep  durations.  The  model  captured  many 
disparate  effects  of  sleep  on  cortisol  dynamics,  such  as  the  inhibition 
of  cortisol  secretion  after  the  wake-to-sleep  transition  and  the  rapid 
rise  of  cortisol  concentration  before  morning  awakening.  Notably,  the 
model  reconciled  the  seemingly  contradictory  hndings  between  stud¬ 
ies  that  report  an  increase  in  cortisol  concentration  following  total 
sleep  deprivation  and  studies  that  report  no  change  in  concentration. 
This  work  provides  a  biomathematical  approach  to  combine  the 
results  on  the  effects  of  sleep  on  cortisol  concentration  into  a  unified 
framework  and  predict  the  impact  of  varying  sleep  durations  on  the 
cortisol  profile. 

biomathematical  models;  hypothalamo-pituitary-adrenal  axis;  sleep 
loss 


CORTISOL  IS  A  KEY  HORMONE  in  the  regulation  of  human  metab¬ 
olism  and  stress  response,  and  its  dysregulation  is  manifested 
in  psychological  disorders  such  as  depression  and  posttrau- 
matic  stress  disorder  (PTSD)  (19,  39)  and  in  metabolic  disor¬ 
ders  such  as  Cushing’s  syndrome  (2).  Cortisol  is  produced  in 
the  zona  fasciculata  of  the  adrenal  cortex.  The  secretion  of 
cortisol  is  regulated  by  adrenocorticotropic  hormone  (ACTH), 
whose  release  from  the  anterior  pituitary  gland  is  induced  by 
the  transport  of  corticotropin-releasing  hormone  (CRH)  from 
the  hypothalamus  to  the  pituitary.  The  secretion  of  ACTH 
occurs  in  pulses  (5,  21),  and  thus  the  secretion  of  cortisol  is 
pulsatile  as  well.  Cortisol  is  distributed  through  the  blood- 
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stream  to  the  hypothalamus  and  pituitary  and  downregulates 
the  secretion  of  CRH  and  ACTH,  thereby  having  a  negative 
feedback  effect  on  its  own  production  (38).  Because  the  secre¬ 
tion  of  CRH  is  governed  by  signals  from  the  suprachiasmatic 
nuclei  in  the  anterior  hypothalamus,  the  basal  cortisol  dynam¬ 
ics  exhibit  a  signihcant  circadian  component  (18,  40).  The  24-h 
cortisol  profile  consists  of  an  early  morning  rise,  decreasing 
levels  during  the  daytime,  and  a  quiescent  period  centered 
around  midnight. 

The  basal  cortisol  dynamics  are  also  affected  by  sleep-wake 
schedules  and  sleep-wake  transitions  (6,  29).  Disturbed  or 
irregular  sleep  schedules  can  dysregulate  cortisol  production, 
and  the  cumulative  effect  of  disturbed  sleep  can  contribute  to 
the  buildup  of  allostatic  load  (25).  For  example,  in  military 
settings,  where  irregular  sleep  schedules  are  common  due  to 
operational  constraints,  irregular  sleep  can  affect  cortisol  levels 
and  have  negative  impacts  on  cognitive  performance,  mood, 
and  stress  (23). 

A  collection  of  mechanistic  models  has  been  proposed  to 
describe  the  various  processes  of  the  hypothalamo-pituitary- 
adrenal  (HPA)  axis,  including  the  neural  firing  of  CRH  (13), 
the  effect  of  glucocorticoid  receptor  count  on  HPA  axis  sta¬ 
bility  (17),  HPA  axis  robustness  to  variations  in  cortisol  bind¬ 
ing  affinity  (20),  and  the  effect  of  variations  in  the  strength  of 
negative  feedback  in  depression  and  PTSD  (31).  However,  the 
existing  mechanistic  models  of  the  HPA  axis  do  not  account 
for  the  effects  of  sleep  on  cortisol  regulation,  which  are  not 
well  understood  at  the  molecular  level  (11).  Instead,  most  of 
the  understanding  of  the  effects  of  sleep  duration  on  cortisol 
production  is  at  a  qualitative,  phenomenological  level  (2).  In 
this  article,  we  present  a  phenomenological  model  describing 
the  effects  of  sleep  duration  on  cortisol  concentration,  thereby 
bringing  together  many  disparate  results  connecting  sleep  and 
cortisol  in  a  unified,  quantitative  framework.  Our  proposed 
model  takes  into  account  the  pulsatile  nature  of  cortisol  secre¬ 
tion,  the  negative  autoregulation  of  cortisol,  and  the  circadian 
and  sleep  allostatic  effects  on  cortisol  concentration. 

We  based  our  model  on  the  Borbely  (3)  two-process  model 
of  sleep  regulation,  where  the  two  processes  are  a  circadian 
process  and  a  sleep  homeostatic  process.  The  sleep  homeo¬ 
static  process  has  its  physiological  basis  in  the  power  generated 
by  electroencephalographic  (EEG)-8  waves  during  slow-wave 
sleep.  As  sleep  progresses,  the  power  of  successive  EEG-8 
wave  episodes  decreases  exponentially.  EEG-8  wave  power  is 
considered  a  measure  of  “sleep  intensity,”  (4)  and  thus,  in  the 
Borbely  (3)  model,  the  sleep  homeostatic  process  is  modeled  as 
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a  decreasing  exponential  during  sleep.  Conversely,  there  is  a 
negative  correlation  between  EEG-8  power  and  the  rate  of 
cortisol  secretion  during  sleep  (16);  thus,  we  hypothesized  that 
the  rate  of  cortisol  secretion  should  increase  exponentially  to  a 
saturation  point  during  sleep.  There  is  also  a  positive  correla¬ 
tion  during  wakefulness  between  EEG-|3  power  and  the  rate  of 
cortisol  secretion  (10).  Because  the  increase  in  EEG-P  power 
corresponds  to  the  effects  of  sleep  deprivation  (24),  we  hy¬ 
pothesized  that  the  rate  of  cortisol  secretion  during  wakeful¬ 
ness  follows  a  saturating  rising  exponential  curve,  as  in  the 
Borbely  (3)  model.  Because  our  hypothesized  process  describ¬ 
ing  the  effects  of  sleep  timing  on  cortisol  secretion  differs  from 
the  Borbely  (3)  sleep  homeostatic  process,  we  refer  to  our 
process  as  the  sleep  allostatic  process. 

We  validated  the  proposed  model  using  data  from  two  studies 
(22,  30)  in  which  one  study  group  experienced  total  sleep  depri¬ 
vation,  a  second  group  experienced  8  h  of  sleep,  and  a  third  group 
experienced  both  sleep  restriction  and  sleep  extension  scenarios. 
We  show  that  our  model  quantitatively  captures  well-established 
relationships  between  sleep  and  cortisol,  such  as  the  inhibition  of 
cortisol  secretion  shortly  after  the  wake-to-sleep  transition  (36) 
and  the  sharp  increase  in  cortisol  concentration  shortly  before 
normal  waking  (29).  Eurthermore,  our  model  reconciles  divergent 
findings  on  the  effects  of  sleep  deprivation  on  cortisol  concentra¬ 
tion  by  demonstrating  when  total  sleep  deprivation  causes  an 
increase  in  cortisol  concentration  (9,  22)  and  when  it  causes  no 
such  increase  (15,  28). 

MATERIALS  AND  METHODS 

This  work  is  a  retrospective  analysis  of  data  originally  reported  by 
Leproult  et  al.  (22)  and  Spiegel  et  al.  (30).  The  pai'ticipants  in  both 
studies  were  healthy  young  men.  Both  protocols  were  approved  by  the 
Institutional  Review  Board  at  the  University  of  Chicago,  and  all 
participants  gave  written  informed  consent. 

Leproult  study  (groups  A  and  B).  Before  the  start  of  the  study,  all 
participants  were  habituated  to  the  laboratory  environment  by  spending 
two  nights  in  the  Clinical  Research  Center  at  the  University  of  Chicago. 
The  participants  were  studied  over  a  32-h  period,  starting  from  1800  on 
day  1  until  0200  on  day  3.  The  participants  were  aware  of  local  clock 
time.  The  participants  remained  recumbent  throughout  the  study  and 
were  maintained  in  dim  light  during  wake  periods  and  in  complete 
darkness  during  sleep  periods.  Food  intake  was  replaced  by  an  intrave¬ 
nous  glucose  infusion  at  a  constant  rate  of  5  g/kg  every  24  h. 

Group  A  consisted  of  17  individuals  [20-30  yr  old,  body  mass 
index  (BMI)  means  ±  SE  22.7  ±  0.5  kg/m^]  who  experienced  total 
sleep  deprivation  during  the  32-h  study  period.  Group  B  consisted  of 
nine  individuals  (22-32  yr  old,  BMI  means  ±  SE  22.8  ±  1.0  kg/m^) 


who  experienced  8  h  of  time  allocated  for  sleep  (TAS)  from  2300  to 
0700.  (Because  the  participants  in  the  2  studies  were  recumbent 
throughout,  we  use  the  phrase  “time  allocated  for  sleep”  instead  of  the 
more  frequently  used  “time  in  bed.”)  Both  groups  were  wakened  at 
0700.  A  sterile  heparin  lock  catheter  was  inserted  in  each  individual’s 
forearm  at  1400,  and  starting  at  1800,  1-ml  blood  samples  were  drawn 
at  20-min  intervals  for  32  h.  The  intravenous  line  was  kept  patent  with 
a  slow  drip  of  heparinized  saline.  Plasma  cortisol  levels  were  deter¬ 
mined  using  the  Coat-A-Count  kit  (Diagnostic  Products,  Los  Angeles, 
CA).  The  lower  limit  of  sensitivity  was  13.8  nmol/1.  The  intra-assay 
coefficient  of  variation  averaged  5%.  All  samples  from  the  same 
individual  were  analyzed  in  the  same  assay.  Eigure  1,  A  and  B,  shows 
a  schematic  of  the  Leproult  study. 

Spiegel  study  (group  C).  Group  C  consisted  of  11  individuals  (18-27 
yr  old,  BMI  means  ±  SE  23.4  ±  0.5  kg/m^).  During  the  week  prior  to  the 
study,  participants  were  asked  to  conform  to  fixed  bedtimes  (2300-0700) 
and  mealtimes.  Wrist  activity  was  monitored  to  verify  compliance.  The 
subjects  spent  16  consecutive  nights  in  the  Clinical  Research  Center, 
consisting  of  three  nights  of  8-h  TAS  (2300-0700),  six  nights  of  4-h 
TAS  (0100-0500),  and  seven  nights  of  12-h  TAS  (2100-0900).  During 
the  last  60  h  of  each  TAS  condition,  the  participants  remained  recumbent. 
During  the  last  24  h  of  the  4-h  TAS  and  12-h  TAS  conditions,  blood  was 
sampled  at  10-  to  30-min  intervals  starting  at  0900.  Participants  received 
identical  carbohydrate-rich  meals  (30  kcal/kg  body  wt,  62%  carbohy¬ 
drates)  at  0900,  1400,  and  1900  during  data  collection.  Plasma  cortisol 
levels  were  measured  by  RIA  (Orion  Diagnostica,  Espoo,  Finland),  with 
a  sensitivity  of  20.7  nmol/1  and  a  4%  intra-assay  coefficient  of  variation. 
Figure  1C  shows  a  schematic  of  the  Spiegel  et  al.  (30)  study.  Subject  8  in 
group  C  required  replacement  of  the  catheter  during  the  4-h  TAS  condition 
and  experienced  a  stress-related  increase  in  cortisol  concentration.  As  a  result 
of  this  confound,  this  subject  was  removed  from  the  analysis. 

Mathematical  analysis.  All  calculations,  parameter  estimations,  and 
cross-validations  were  performed  in  MATLAB  R201  IB.  Nonlinear  least 
squares  estimation  was  used  for  parameter  estimation.  For  the  0-  and  8-h 
TAS  conditions,  pulsatile  cortisol  secretion  was  set  to  begin  at  0700  on 
day  7;  cortisol  secreted  before  that  time  was  assumed  to  have  disappeared 
by  the  time  data  collection  started.  For  the  4-  and  12-h  TAS  scenarios, 
secretion  was  set  to  begin  at  1900  on  days  7  and  14,  respectively. 

The  nonlinear  coefficient  of  determination  (r^)  between  a  fit  f(0 
and  a  data  set  y(ri),  y(f2),  .  .  .  y(fn)  was  calculated  as 

=  1  -  [(2r^i  -  y(fi))")/(2r^i  (y(fi)  -  y)")]. 

where  y  is  the  mean  of  the  data  set.  The  root  mean  squared  error 
(RMSE)  between  the  fit  and  the  data  set  was  calculated  as 

RMSE  =  V(2i=i(f(fi)-y(fi))")/n- 

RESULTS 

Two-process  model  for  cortisol  secretion  and  concentration. 
We  describe  the  rate  of  cortisol  secretion  using  a  two-process 
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Fig.  1.  Schematics  of  the  study  protocols.  Black  bars  indicate  time  allocated  to  sleep  (TAS).  Gray  bars  indicate  times  when  data  collection  occurred.  A:  protocol 
experienced  by  group  A.  B\  protocol  experienced  by  group  B.  C:  protocol  experienced  by  group  C.  During  the  sleep  restriction  nights,  the  TAS  was  from  0100 
to  0500.  During  the  sleep  extension  nights,  the  TAS  was  from  2100  to  0900.  Data  were  collected  for  24  h  starting  at  0900  on  the  8th  and  15th  days,  d-hh,  day-hour 
(indicating  that  the  time  given  on  the  jc-axis  is  given  by  the  day  within  the  study  followed  by  the  2  digits  indicating  the  time  of  day;  e.g.,  1-18  indicates  that 
the  time  is  1800  on  day  1  of  the  study). 
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model  that  contains  both  a  circadian  component  and  a  sleep 
allostatic  component.  The  circadian  function  C(f),  measured  in 
nmol/h,  is  defined  at  time  t  as 

C(r)  =  ajsin 

where  the  ai  denotes  the  Borbely-Achermann  parameters  ai  = 
0.97,  &2  =  0.22,  a.2  =  0.07,  and  a4  =  0.03  (4),  p  represents  the 
circadian  amplitude  in  nmol/h,  and  0  is  the  circadian  phase  in 
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Fig.  2.  Development  of  the  two-process  model  (TPM)  for  cortisol.  A:  the 
circadian  process  C  describes  an  entrained  circadian  rhythm  governed  by 
amplitude  (3  and  phase  6.  B\  the  sleep  allostatic  process  S  describes  the  effect 
of  sleep  on  the  rate  of  cortisol  production.  During  wake,  process  5  is  a  rising 
saturating  exponential  governed  by  amplitude  aw  and  rate  parameter 
During  sleep,  process  S  is  governed  by  amplitude  as  and  rate  parameter  ys. 
C:  the  TPM  generates  pulses  whose  amplitude  is  defined  in  part  by  the  sum  of 
processes  S  and  C.  The  magnitude  of  these  pulses  is  decreased  according  to  a 
proportional  feedback  constant,  ^p.  Each  physiological  pulse  is  modeled  as  a 
pair  of  Kronecker  8-functions  occumng  with  a  period  of  80  min.  D:  plasma 
cortisol  concentration  is  modeled  as  the  sum  of  decaying  exponential  functions 
with  cortisol  disappearance  rate  dc.  The  parameter  values  used  in  this  plot  are 
the  mean  pai'ameter  values  for  group  B,  which  are  shown  in  Table  2. 
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Fig.  3.  A:  group  mean  raw  data  and  model  fit  for  group  A{n  =  17).  B\  group  mean 
raw  data  and  model  fit  for  group  B  {n  =  9).  The  black  bar  indicates  the  8-h  TAS 
period.  Both  groups  slept  from  2300  the  previous  night  to  0700  on  study  day  I.  In 
the  fitted  mean  cortisol  concentration,  the  Mann- Whitney  t/-test  showed  that 
group  B  has  a  significantly  lower  concentration  of  cortisol  from  2000  on  day  2  to 
0000  on  day  3  (P  =  0.006),  corroborating  the  result  in  Ref.  22. 
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Fig.  4.  A:  group  mean  raw  data  and  model  fit  for  group  C  (n  =  10)  in  the  4-h 
TAS  scenario.  Black  bai‘  indicates  the  TAS  period.  B:  group  mean  raw  data  and 
model  fit  for  group  C  in  the  12-h  TAS  scenario.  Black  bar  indicates  the  TAS 
period.  There  was  a  large  spike  in  the  4-h  TAS  data  at  0530  on  day  9  that  was 
not  captured  by  the  model  fit.  We  believe  this  spike  was  due  to  the  cortisol 
awakening  response  (see  discussion).  There  was  also  a  saturation  in  the  12-h 
TAS  data  at  0730  on  day  16  that  was  not  well  captured  by  the  model  fit.  This 
leveling  off  is  likely  due  to  the  majority  of  individuals  waking  while  remaining 
in  sleep  position. 
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Fig.  5.  Representative  individual  fits  for  the  0-  (A),  8-  (B),  4-  (C),  and  12-h  TAS  (D)  scenarios.  The  individual  fits  demonstrated  pulsatile  behavior  in  cortisol 
concentration,  whereas  the  group  average  fits  did  not.  The  individual  fits  chosen  for  this  plot  were  those  with  the  median  root  mean  squared  error  among  the 
fits  for  each  scenario.  Figures  7-10  show  the  full  set  of  individual  fits. 


hours.  The  sleep  allostatic  function  S(t),  also  measured  in 
nmol/h,  is  piecewise  continuous  and  consists  of  two  saturating 
rising  exponentials: 

“w(l  “  during  wake, 

“s(l  “  Is  during  sleep, 

where  the  parameters  during  wake  are  Tsw,  the  most  recent 
sleep-to-wake  transition  time,  ct*,  the  wake  allostatic  magni¬ 
tude  in  nmol/h,  and  -y^,  the  unitless  wake  allostatic  rate.  During 
sleep,  the  parameters  are  T^s,  the  most  recent  wake-to-sleep 
transition  time,  ots,  the  sleep  allostatic  magnitude  in  nmol/h, 
and  js,  the  unitless  sleep  allostatic  rate.  S(t)  models  a  rate  of 
secretion,  which  changes  discontinuously  at  each  sleep-wake 
transition  when  its  value  instantaneously  decreases  to  zero.  The 
two-process  model  output  at  time  t  is 
00 

TPM(T)  =  J  [S(t)  -f  C(?)]8(?  -  T)dr, 

—  00 

where  8  is  a  Dirac  8-function.  TPM(T)  is  measured  in  nmol/1. 

Cortisol  is  produced  in  pulses.  We  assumed  an  interpulse 
period  of  80  min  (13)  and  modeled  each  physiological  pulse  as  a 
pair  of  Kronecker  8-functions  occurring  10  min  apart.  This  al¬ 
lowed  us  to  capture  the  situations  where  a  cortisol  measurement 
occurs  before,  after,  or  during  a  20-min-long  pulse.  We  assumed 


that  the  phase  of  the  pulsatile  rhythm  is  locked  to  the  circadian 

1 

phase,  with  8-functions  occurring  at  times  0  -I-  2n  4-  —  and  0  -I- 
3 

2n  +  — ,  where  n  is  any  integer  number  of  hours.  We  denoted  the 

set  of  times  at  which  8-functions  occur  as  Tg. 

The  concentration  of  plasma  cortisol  y(t)  at  time  t  is  defined  as 

y(0  =  S  [TPM(T)  -  kpy(T)]e-‘'c('-T)^ 

T<r,TeT8 

where  kp  is  a  proportional  feedback  constant  describing  the 
effect  of  cortisol  autoregulation  and  dc  is  the  cortisol  disap¬ 
pearance  rate.  The  concentration  y{t)  is  modulated  primarily  by 
the  changing  amplitude  of  the  secretion  pulses  (35).  We  thus 
modeled  an  individual’s  plasma  cortisol  concentration  as  a  sum 
of  decaying  exponentials  with  eight  parameters:  circadian  am¬ 
plitude  (3  and  phase  0,  wake  allostatic  amplitude  otw  and  rate 
7w,  sleep  allostatic  amplitude  as  and  rate  js,  cortisol  feedback 
constant  kp,  and  disappearance  rate  dc.  The  constituent  ele¬ 
ments  of  the  model  are  illustrated  in  Fig.  2. 

To  evaluate  the  performance  of  the  model,  we  first  fit  the 
model  to  individual  and  group  mean  data  for  different  sleep 
timing  scenarios.  We  then  determined  whether  the  parameter 
values  calculated  for  individuals  in  each  study  group  are 


Table  1.  Group  mean  and  individual  measurements  of  goodness  of  fit  for  each  of  the  4  study  scenarios 


TAS,  h 

Group  Mean  RMSE 

Mean  (SD)  of  Individual  Fit  RMSE 

Group  Mean  /•^ 

Mean  (SD)  of  Individual  Fit 

0 

34.2 

80.3  (21.9) 

0.91 

0.66(0.10) 

4 

46.1 

67.6  (16.7) 

0.79 

0.64  (0.09) 

8 

36.7 

78.6  (20.8) 

0.89 

0.65  (0.11) 

12 

37.4 

67.6(17.5) 

0.91 

0.75  (0.12) 

TAS,  time  allocated  to  sleep;  RMSE,  root  mean  squared  enor  (nmol/1);  r^,  coefficient  of  determination. 
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Table  2.  Means  (SD)  of  the  fitted  individualized  model  parameters  for  each  of  the  3  study  groups 


Study  Group  aw,  nmol/h  (3,  nmol/h  0  (Time)  dc,  1/h  A:p  as,  nmol/h  ys 


Group  A  70.8  (34.6)  0.954  (0.025)  87.4  (17.2)  0239  (0042)  0.448  (0.061)  0.246  (0.060) 

Group  B  44.2  (24.2)  0.917  (0.034)  64.1  (19.0)  0233  (0017)  0.422  (0.061)  0.169  (0.123)  173.2  (79.1)  0.895  (0.035) 

Group  C  26.0(9.18)  0.902(0.036)  46.4(8.58)  0321  (0026)  0.331  (0.031)  0.148  (0.070)  103.0(51.6)  0.951  (0.028) 


ttw,  sleep-to-wake  transition  time;  unitless  wake  allostatic  magnitude  in  nmol/h;  p,  circadian  amplitude;  6,  circadian  phase;  dc,  cortisol  disappearance  rate; 
kp,  cortisol  feedback  constant;  as,  sleep  allostatic  magnitude  in  nmol/h;  ys,  unitless  sleep  allostatic  rate. 


consistently  distributed  across  sleep  scenarios.  Finally,  we 
cross-validated  across  study  scenarios  by  comparing  the  fit 
generated  by  one  group  in  a  given  scenario  to  the  fit  of  the 
cortisol  profile  predicted  by  our  model  for  a  different  study 
group  undergoing  the  same  scenario. 

Individual  and  group  mean  model  fits.  For  each  of  the  three 
study  groups,  we  obtained  group  mean  model  fits  by  first 
calculating  a  fit  for  each  individual  in  the  group  and  then 
averaging  the  individualized  fits  to  determine  the  group  mean 
fit.  We  used  this  procedure  because  the  pulsatile  secretions  of 
cortisol  are  readily  apparent  in  the  individual  data  but  obscured 
in  the  group  mean  data  because  of  between-subject  variations 
in  the  phase  6.  For  group  C,  we  calculated  its  individuals’  fits 
by  simultaneously  fitting  the  data  from  both  scenarios  to 
generate  one  set  of  parameter  values  for  each  individual. 

Figure  3  shows  the  group  mean  fits  for  group  A  in  the  0-h  TAS 
scenario  and  group  B  in  the  8-h  TAS  scenario.  In  the  fitted  mean 
cortisol  concentration,  the  Mann- Whitney  [/-test  showed  that 
group  B  has  a  significantly  lower  concentration  of  cortisol  from 
2000  on  day  2  to  0000  on  day  3  (P  =  0.006),  corroborating  the 
result  in  Ref.  22.  Figure  4  shows  the  group  mean  fits  for  group  C 
in  the  4-  and  12-h  TAS  scenarios.  Figure  5  shows  representative 
individualized  fits  for  each  of  the  four  scenarios,  which  illustrate 
the  pulsatile  nature  of  cortisol  concentration. 

Table  1  shows  the  goodness  of  fit  for  the  individualized  fits 
and  for  each  group  mean  fit,  measured  in  terms  of  RMSE  and 
r^.  The  individual  RMSE  values  were  smaller  for  the  4-  and 
12-h  TAS  scenarios,  whereas  the  group  RMSE  values  were 
smaller  for  the  0-  and  8-h  TAS  scenarios.  Averaging  over  the 
four  scenarios,  the  r^  values  indicated  that  the  model  accounted 
for  88%  of  the  variance  in  the  group  mean  data  and  67%  of  the 
variance  in  the  individual  data. 

Comparison  of  model  parameters  between  groups.  Table  2 
shows  the  mean  and  standard  deviation  of  each  model  param¬ 
eter  in  each  study  group.  The  sleep  allostatic  parameters  as  and 
js  could  not  be  determined  for  group  A  because  they  did  not 
sleep  in  their  scenario.  Table  3  shows  the  P  values  calculated 
from  Kolmogorov-Smirnov  tests  on  the  distributions  of  indi¬ 
vidualized  parameters  between  study  groups.  Excluding  the 
phase  parameter  0,  which  is  a  state  parameter  (27)  that  depends 
on  sleep-wake  history  and  other  environmental  conditions,  there 
are  no  statistically  significant  differences  between  the  group  B 
parameters  and  those  of  groups  A  and  C.  However,  the  parameters 
for  groups  A  and  C  are  significantly  different.  The  significant 
differences  between  groups  A  and  C  suggest  that  the  model  is 
insensitive  to  the  absolute  amplitudes  of  the  model’s  produc¬ 
tion  parameters  (aw,  (3,  as)  and  degradation  parameters  (kp,  dc) 
because  the  estimates  for  all  of  these  parameters  are  higher  for 
group  A  than  for  group  C.  The  model  is  more  sensitive  to  the 
ratios  between  production  and  feedback  parameters.  Between 
groups  A  and  C,  P  values  calculated  from  Kolmogorov- 


Smirnov  tests  indicate  that  the  distributions  of  the  ratios 
a^/kp  (P  =  0.18)  and  (BZ/cp  (P  =  0.25)  are  not  significantly 
different  between  the  groups. 

Cross-validation  between  study  scenarios.  We  performed 
cross-validation  by  substituting  the  parameters  for  the  individ¬ 
uals  in  group  C  into  the  models  for  the  0-  and  8-h  TAS 
scenarios  and  by  substituting  the  parameters  for  the  individuals 
in  group  B  into  the  models  for  the  0-,  4-,  and  12-h  TAS 
scenarios.  We  did  not  perform  cross-validation  from  group  A 
onto  the  other  study  scenarios  because  sleep  allostatic  param¬ 
eters  were  not  available.  To  account  for  differences  in  cortisol 
amplitude  across  studies,  we  multiplied  cortisol  concentrations 
by  1.14  when  cross-validating  group  C  on  the  0-  and  8-h  TAS 
scenarios  and  divided  by  1.14  when  cross-validating  group  B 
on  the  4-  and  12-h  TAS  scenarios.  We  determined  the  value 
1.14  by  taking  the  ratio  of  the  grand  mean  cortisol  level  of 
groups  A  and  B  combined  with  that  of  group  C.  To  account  for 
differences  in  phase,  we  shifted  the  phase  parameter  6  by 
—  1.02  h  when  cross-validating  group  B  on  the  0-h  TAS 
scenario,  by  —0.77  h  when  cross-validating  group  B  on  the  4- 
and  12-h  TAS  scenarios,  by  —0.25  h  when  cross-validating 
group  C  on  the  0-h  TAS  scenario,  and  by  0.77  h  when 
cross-validating  group  C  on  the  8-h  TAS  scenario. 

Figure  6  shows  the  adjusted  group  mean  plasma  cortisol  fits. 
For  model  parameters  from  groups  B  and  C,  the  Mann- 
Whitney  [/-test  showed  that  the  value  of  the  fit  onto  the  8-h 
TAS  scenario  is  significantly  greater  than  the  value  of  the  fit 
onto  the  0-h  TAS  scenario  from  2000  on  day  2  to  0000  on  day 


Table  3.  P  values  resulting  from  2-sample 
Kolmogorov-Smirnov  tests  on  the  distributions 
of  parameter  values  calculated  for  each  study  group 


Model  Parameters 

Comparisons  of  Study  Groups 

Group  A  vs.  B 

Group  A  vs.  C 

Group  B  vs.  C 

aw 

0.21 

0.00012 

0.07 

7w 

0.017 

0.0036 

0.45 

P 

0.017 

0.0000064 

0.015 

0 

0.00022 

0.18 

0.0026 

dc 

0.58 

0.0000290 

0.0091 

kp 

0.073 

0.00078 

0.9 

as 

0.044 

Js 

0.0091 

Significant  differences  are  in  boldface.  The  Bonferroni-corrected  threshold  of 
significance  is  P  =  0.05  -i-  8.  The  significant  differences  in  the  distribution  of  6  can 
be  attributed  to  adjustments  in  the  circadian  rhythm  following  6  days  of  sleep 
restriction/extension  or  possibly  seasonal  differences  in  phase.  The  significant 
differences  in  the  distributions  of  the  parameters  for  groups  A  and  C  suggest  that 
the  model  is  insensitive  to  the  absolute  amplitudes  of  the  model’s  production  (aw, 
(3,  and  as)  and  degradation  parameters  {kp,  dc),  because  the  estimates  for  all  of 
these  parameters  are  higher  for  group  A  than  for  group  C. 
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Fig.  6.  Adjusted  cross-validation  tits  between  the  different  study  scenarios:  0-  (A),  8-  (B),  4-  (C),  and  12-h  TAS  (D).  The  fits  were  adjusted  with  an  amplitude 
adjustment  factor  of  1.14  to  account  for  differences  in  study  conditions  and  phase  adjustments  to  account  for  differences  in  the  group  mean  estimates  of  phase. 
Black  line  in  A  indicates  the  fit  calculated  from  the  group  A  parameters,  the  dark  gray  lines  in  A-D  the  fit  from  the  group  B  parameters,  and  the  light  gray  lines 
in  in  A-D  the  fit  from  the  group  C  pai'ameters.  Solid  lines  indicate  self-fits,  and  dashed  lines  indicate  cross-validation  fits.  Black  bars  indicate  the  TAS.  Figure 
1 1  shows  the  unadjusted  cross-validations.  For  model  parameters  from  groups  B  and  C,  the  Mann-Whitney  f/-test  showed  that  the  value  of  the  fit  onto  the  8-h 
TAS  scenario  is  significantly  greater  than  the  value  of  the  fit  onto  the  0-h  TAS  scenario  from  2000  on  day  2  to  0000  on  day  3  (P  =  0.003  for  group  B,  P  = 
0.0001  for  group  C),  in  accord  with  the  within-subject  data  reported  in  Ref.  9. 


3  (P  =  0.003  for  group  B,  P  =  0.0001  for  group  C),  in  accord 
with  the  within-suhject  data  reported  in  Ref.  9. 

Table  4  shows  the  goodness  of  the  cross-validation  fits  in 
terms  of  r^,  RMSE,  and  the  P  value  from  the  Mann-Whitney 
[/-test,  which  measures  the  probability  that  the  ht  has  the  same 
median  as  the  data.  The  r^  values  for  the  cross-validation  hts 
ranged  from  52  to  90%,  and  the  RMSE  values  ranged  from 
35.5  to  85.5  nmol/1.  The  relatively  poor  goodness  of  fit  for  the 
cross-validation  from  group  B  onto  the  12-h  TAS  scenario  was 
due  to  the  fact  that  individuals  in  group  C  did  not  sleep  for  the 
whole  12  h;  the  average  time  asleep  for  individuals  in  the 
scenario  was  9  h  and  3  min  (30). 

DISCUSSION 

We  presented  a  phenomenological  model  that  describes  the 
effects  of  sleep  duration  on  plasma  cortisol  concentration.  We 
based  the  model  on  the  Borbely  two-process  model  of  sleep 


regulation  (3)  and  defined  the  amplitude  of  cortisol  pulses  as 
the  combination  of  a  sleep  allostatic  process  S,  a  circadian 
process  C,  and  a  negative  feedback  term  kp.  The  structure  of 
process  S  was  inferred  from  the  results  of  Gronfier  et  al.  (16), 
which  show  a  negative  correlation  between  cortisol  secretion 
rate  and  EEG-(3  power  during  sleep,  and  the  results  of  Chapo- 
tot  et  al.  (10),  which  show  a  positive  correlation  between  the 
secretion  rate  and  EEG-|3  power  during  wake.  Process  C  was 
constructed  using  the  standard  parameters  used  Borbely  and 
Achermann  (4)  for  modeling  sleep  regulation.  The  proportional 
feedback  term  kp  is  a  simplified  representation  of  cortisol’s 
negative  feedback  mechanism  (38)  in  which  cortisol  down- 
regulates  the  production  of  CRH  and  ACTH. 

The  phenomenological  model  predicts  the  decrease  of  cortisol 
concentration  observed  after  the  wake-to-sleep  transition  by 
Weitzman  et  al.  (36)  and  the  rapid  increase  in  cortisol  concentra¬ 
tion  before  normal  waking  observed  by  Spath-Schwalbe  et  al. 


Table  4.  Measurements  of  goodness  of  fit  for  the  cross-validations 


TAS,  h 


Itudy  Group 

0 

4 

8 

12 

1^ 

RMSE 

P 

RMSE 

P 

1^ 

RMSE 

P 

jO. 

RMSE 

P 

Group  A 

0.91 

34.2 

0.67 

Group  B 

0.90 

35.5 

0.78 

0.74 

50.2 

0.96 

0.89 

36.7 

0.90 

0.52 

85.5 

0.04 

Group  C 

0.88 

39.3 

0.54 

0.78 

45.7 

0.90 

0.68 

61.9 

0.58 

0.92 

35.3 

0.73 

Before  these  measurements  were  calculated,  the  fits  were  adjusted  with  an  amplitude  adjustment  factor  of  1.14  and  phase  adjustments  to  account  for  differences 
in  the  group  mean  estimates  of  phase.  No  adjustments  were  made  to  the  cross-validation  fits  within  studies.  Values  in  boldface  indicate  self-fits.  No 
cross-validations  were  made  from  group  A  onto  other  scenarios  because  the  sleep  allostatic  parameters  were  unavailable.  P  is  the  P  value  from  Mann-Whitney 
f/-test.  BonfeiToni-corrected  threshold  for  significance:  P  <  0.05  ^  9.  Table  5  shows  the  values  of  these  measurements  for  the  unadjusted  cross-validations. 
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Fig.  7.  Fits  for  each  of  the  individuals  in  group  A  in  the  0-h  TAS  scenario.  Subject  1  had  the  median  root  mean  squared  error  among  the  17  individuals,  and 
the  plot  for  subject  1  from  Fig.  7  is  identical  to  the  plot  in  Fig.  5A. 


(29).  Both  of  these  phenomena  are  explained  by  the  exponential 
form  of  Process  S  during  sleep.  At  each  wake-to-sleep  transition, 
the  sleep  allostatic  contribution  to  cortisol  secretion  drops  to  zero, 
leading  to  the  decrease  in  plasma  cortisol  for  1-2  h  after  sleep 
onset  observed  by  Weitzman  et  al.  (36).  Also,  the  estimated  values 
of  the  sleep  allostatic  magnitude  as  predicted  a  rapid  increase  in 
the  allostatic  contribution  to  cortisol  secretion  after  5-8  h  of  sleep. 


leading  to  the  rapid  increase  in  concentration  observed  by  Spath- 
Schwalbe  et  al.  (29). 

Notably,  our  model  reconciles  reports  that  there  is  no  change  in 
cortisol  concentration  as  a  result  of  total  sleep  deprivation  (15, 
28),  with  reports  claiming  that  there  is  an  increase  (9,  22).  The 
cause  of  the  apparent  discrepancy  is  the  time  of  measurement. 
Follenius  et  al.  (15)  and  Salm-Pascual  et  al.  (28)  measure  the  night 
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Fig.  8.  Fits  for  each  of  the  individuals  in  group  B  in  the  8-h  TAS  scenario.  Subject  1  had  the  median  root  mean  squared  en'or  among  the  9  individuals,  and  the 
plot  for  subject  1  from  Fig.  7  is  identical  to  the  plot  in  Fig.  5B.  Black  bars  indicate  the  TAS. 


cortisol  profile  from  2300  to  0700  and  from  2200  to  0600, 
respectively.  During  these  hours,  our  model  predicted  that  8-h 
TAS  individuals  would  experience  a  small  decrease  in  cortisol 
concentration  immediately  after  sleep  onset  as  the  contrihution  of 
the  allostatic  process  S  decreased  to  zero.  However,  because  the 
sleep  allostatic  amplitude  as  was  greater  than  the  wake  allostatic 
amplitude  a^  (Table  2),  the  cortisol  concentration  of  the  sleeping 
individuals  increased  more  rapidly  than  those  of  the  awake  indi¬ 
viduals,  leading  to  no  significant  total  difference  over  the  8-h  sleep 
period.  In  fact,  our  model  predicted  a  slightly  lower  cortisol 
concentration  in  8-h  TAS  individuals  just  after  sleep  onset  and  a 
slightly  higher  concentration  just  before  waking,  qualitatively 
mimicking  the  data  of  Salm-Pascual  et  al.  (28).  Whereas  Leproult 
et  al.  (22)  compared  a  0-h  TAS  group  and  an  8-h  TAS  group  and 
Chapotot  et  al.  (9)  reported  on  a  within-subject  study  on  the  same 
group,  both  studies  measured  cortisol  concentration  for  at  least  1 1 
h  after  waking.  Our  model  predicted  that  differences  in  cortisol 
between  groups  A  and  B  occur  in  the  evening  after  sleep  depri¬ 
vation  (Fig.  3).  Also,  we  cross-validated  group  B  onto  the  0-h 
TAS  scenario  and  observed  that  the  predicted  cortisol  levels  the 
following  evening  are  greater  than  in  the  normal  sleep  scenario,  in 
agreement  with  the  within-subject  results  of  Chapotot  et  al.  (9) 
(Fig.  6). 


To  minimize  the  number  of  model  parameters,  we  made 
several  simplifying  assumptions.  First,  we  assumed  a  phase 
lock  between  the  phase  of  cortisol  secretion  pulses  and  the 
phase  of  the  circadian  process  C.  Second,  we  assumed  a 
simplified  version  of  the  cortisol  autoregulation  process.  In¬ 
stead  of  explicitly  modeling  the  feedback  mechanisms  by 
which  cortisol  represses  its  precursors  CRH  and  ACTH  (38), 
we  modeled  the  negative  feedback  loop  with  a  single  propor¬ 
tional  constant,  kp.  Third,  we  assumed  a  fixed  interpulse  period 
of  80  min,  although  this  period  can  vary  within  and  between 
individuals.  By  removing  or  modifying  these  assumptions,  we 
could  produce  a  more  comprehensive  but  less  parsimonious 
model. 

The  two-process  model  parameters  show  significant  be- 
tween-subject  variability  (Table  2).  Some  of  this  variability 
is  due  to  the  difficulty  in  estimating  absolute  values  of  both 
production  and  degradation  parameters  (Table  3),  but  there 
is  significant  variability  between  subjects  that  is  caused  by 
underlying  physiological  factors.  Further  investigation  is 
necessary  to  relate  the  parameters  from  our  model  to  the 
neuroendocrine  parameters  from  mechanistic  models  (5,  13, 
17,  20,  31).  We  hypothesize  that  between-subject  variability 
in  glucocorticoid  receptor  counts  (34)  could  account  for  the 
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Fig.  9.  Fits  for  each  of  the  individuals  in  group  C  in  the  4-h  TAS  scenario.  Subject  2  had  the  6th  largest  root  mean  squared  error  among  the  10  individuals,  and 
the  plot  for  subject  2  from  Fig.  7  is  identical  to  the  plot  in  Fig.  5C.  Black  bars  indicate  the  TAS. 


variability  in  the  cortisol  feedback  constant  kp  and  disap¬ 
pearance  rate  dc  and  that  between-subject  variability  in  the 
strength  of  the  pulse  generator  in  the  pituitary  (21)  may  be 
the  cause  of  variability  in  the  amplitude  parameters  a*,  as, 
and  (3.  Differences  in  the  CRH  and  ACTH  secretion  systems 
in  the  hypothalamus  and  pituitary  may  account  for  the 
variability  in  the  rate  parameters  y/w  and  y/s-  A  possible 
confound  affecting  the  distribution  of  the  disappearance  rate 
is  the  posture  of  the  individual;  the  mean  value  we  report  for 
dc  yields  a  mean  cortisol  half-life  of  105  min  (SD  =  19 
min),  slightly  larger  than  the  value  reported  in  (14).  How¬ 
ever,  this  rate  is  likely  affected  by  the  recumbent  posture  of 
the  individuals  in  the  studies  we  considered  (33). 


Furthermore,  the  effects  of  the  cortisol  awakening  response  (26, 
37)  and  of  meals  (33)  are  key  phenomena  not  included  in  our  model 
(1).  The  effect  of  the  cortisol  awakening  response  can  be  seen  as  a 
spike  in  both  group  mean  and  individual  cortisol  concentrations  at  the 
sleep-to-wake  transition  in  the  4-h  TAS  scenario  (Fig.  4),  whereas  the 
effect  of  meals  can  be  seen  in  the  postprandial  increases  in  cortisol  in 
group  C  (Fig.  3).  Despite  not  exphcitly  modeling  these  phenomena 
that  were  observed  only  in  the  4-  and  12-h  TAS  scenarios,  the 
goodness  of  fit  of  the  model  in  terms  of  RMSF  and  was  similar  for 
all  groups.  However,  the  cortisol  awakening  response  in  particular 
has  significant  impact  on  the  diurnal  cortisol  profile,  and  modeling 
this  phenomenon  will  be  essential  to  understanding  how  varying 
sleep  durations  impact  the  time  course  of  cortisol  concentration. 
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Fig.  10.  Fits  for  each  of  the  individuals  in  group  C  in  the  12-h  TAS  scenario.  Subject  2  had  the  6th  largest  root  mean  squared  error  among  the  10  individuals, 
and  the  plot  for  subject  2  from  Fig.  7  is  identical  to  the  plot  in  Fig.  5D.  Black  bars  indicate  the  TAS. 


Further  research  is  also  needed  to  model  differences  in 
cortisol  profiles  due  to  sex  and  age  (32),  sleep  shifting  (7,  8), 
psychological  disorders  such  as  PTSD  and  depression,  and 
metabolic  disorders  such  as  Cushing’s  syndrome  (25).  We 
expect  that  our  phenomenological  model  will  provide  insight 
into  these  causes  of  between-subject  variation  in  cortisol  con¬ 
centration  and  also  have  potential  application  to  many  other 
metabolites  (12)  and  hormones  that  also  display  circadian 
rhythmicity. 

APPENDIX 

Figures  7,  8,  9,  and  10  show  the  individualized  fits  for  cortisol 
concentration  for  the  0-,  8-,  4-,  and  12-h  TAS,  respectively. 

Figure  11  shows  the  unadjusted  group  mean  plasma  cortisol  fits 
wherein  we  did  not  account  for  differences  in  cortisol  amplitude  and 
circadian  phase.  Table  5  shows  the  goodness  of  the  cross-validation 


fits  in  terms  of  r^,  RMSE,  and  the  P  value  from  the  Mann- Whitney 
[/-test.  The  values  for  the  cross-validation  fits  ranged  from  —  1  to 
83%,  and  the  RMSE  values  ranged  from  46.8  to  123.5  nmol/1.  The 
poor  goodness  of  fit  for  the  cross-validation  from  group  B  onto  the 
12-h  TAS  scenario  was  due  to  the  fact  that  individuals  in  group  C  did 
not  sleep  for  the  whole  12  h,  the  difference  in  circadian  phase  between 
group  B  and  group  C,  and  the  greater  mean  cortisol  amplitude  of 
groups  A  and  B  compared  with  group  C. 
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Fig.  1 1 .  Unadjusted  cross-validation  fits  between  the  different  study  scenarios:  0-  (A),  8-  (B),  4-  (C),  and  12-h  TAS  (D).  Black  line  in  A  indicates  the  fit  calculated 
from  the  group  A  parameters,  the  dark  gray  lines  in  A-D  the  fit  from  the  group  B  parameters,  and  the  light  gray  lines  in  A-D  the  fit  from  the  group  C  parameters. 
Solid  lines  indicate  self-fits,  and  dashed  lines  indicate  cross-validation  fits.  Black  bars  indicate  the  TAS  for  each  scenario.  Figure  6  shows  the  adjusted 
cross-validations. 
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